#!/usr/bin/env python
# -*- coding: utf-8 -*-
from moog import moog
import pylab, atpy,param_finder_modules
import numpy as np
import os

def spectrum_loader(nS,specFile):
	words=specFile.split()
	ifile = open('/home/jeffrey/PHD/ruawahia/data/jaccodata/'+words[nS-1],'r')
	wave=[]; intensity=[]
	for line in ifile:
	    quadruple = line.split()
	    if quadruple[1]=="*********":
			continue
	    wave.append(float(quadruple[0])); intensity.append(float(quadruple[2]))
	return wave, intensity, words[nS-1]

params_found = open('params_news_out','w')
tbl = atpy.Table('data/v07_WFI_2MASS.vot')
for i in range(len(tbl)):
	#tbl = tbl_overall.where((tbl_overall.LEID == LEID))
	LEID=tbl.data['LEID'][i]
	print LEID
	RA=tbl.data['RAJ2000'][i]
	Dec=tbl.data['DEJ2000'][i]
	Vmag=tbl.data['Vmag_B09'][i]
	Bmag=tbl.data['Bmag_B09'][i]
	nS=tbl.data['nS'][i]
	spFile=tbl.data['spFile'][i]
	J_2MASS=tbl.data['j_m'][i]
	K_2MASS=tbl.data['k_m'][i]
	FeH_initial=-1.4999
	BmV=(Bmag-Vmag)-0.12
	Vmag=Vmag-0.36 #This is the A(V)		
	temp_BmV=param_finder_modules.Alonso_BmV(BmV,FeH_initial)
	
	#The A(V) is unapplied since the 2MASS photometry has not been corrected
	[JmK,VmK]=param_finder_modules.TCS_transforms(J_2MASS,K_2MASS,Vmag+0.36)
	print JmK,VmK
	temp_VmK=param_finder_modules.Alonso_VmK(VmK,FeH_initial)
	M_V=Vmag-13.7
	if temp_VmK=="":
		continue
	log_g=param_finder_modules.Alonso_logg(temp_VmK,FeH_initial,M_V)
	print temp_VmK, log_g
	temp_round=int(round(temp_VmK/50)*50)
	grav_round=round(log_g/0.1)*0.1
	
	def find_nearest(array,value):
	    idx=(np.abs(array-value)).argmin()
	    return idx
	
	#Load the vL07 spec and convert to a numpy array
	[wave,intensity,holder]=spectrum_loader(1,spFile)
	jacco_spec = np.array([wave,intensity])
	jacco_spec = np.transpose(jacco_spec)
	
	chi2_best=99
	for FeH in range(24):
		FeH =float(FeH)/-10-0.2
		#Use MOOG to compute the spectrum
		model_atm="moog_holder.model"
		linelists="/home/jeffrey/PHD/software/MOOG2010/linelists/lin4500"
		for input_ba in range(13):
			input_ba=float(input_ba)/10-0.2
			spectrum=[]
			try:
				spectrum=np.load("/media/store/MOOG_synspec/"+str(temp_round)+"_"+str(grav_round)+"_"+str(FeH)+"_"+str(input_ba)+".npy")
			except (IOError):
				#Create the moog model using atmosphy
				os.system("atmosphy -teff "+str(temp_round)+" -logg "+str(grav_round)+" -feh "+str(FeH)+" -o moog_holder.model")
				spectrum=moog.run(model_atm,linelists,None,[4500,4600],ba_abun=input_ba)
				os.remove("moog_holder.model")
				np.save("/media/store/MOOG_synspec/"+str(temp_round)+"_"+str(grav_round)+"_"+str(FeH)+"_"+str(input_ba), spectrum)
			#This section matches the two spectra together and produces two arrays containing just the matches
			#It has to round the vL07 spectra as the wavelength points are more finely given than MOOG calculates
			c = np.in1d(np.round(jacco_spec[:,0],2),spectrum[:,0])
			jacco_match=jacco_spec[c,:]
			
			c = np.in1d(spectrum[:,0],np.round(jacco_spec[:,0],2))
			spectrum_match=spectrum[c,:]
			
			
			cont_points=np.array([4538.00,4559.50,4567.50])
			cont_diff=float()
			for i in range(3):
				cont_idx=find_nearest(jacco_match[:,0],cont_points[i])
				#cont_idx_moog=find_nearest(jacco_spec[:,i],cont_points[i])
				cont_diff=jacco_match[cont_idx,1]-spectrum_match[cont_idx,1]+cont_diff
			cont_shift=cont_diff/3
			
			spectrum_match[:,1]=spectrum_match[:,1]+cont_shift
			#pylab.plot(spectrum_match[:,0],spectrum_match[:,1],'-')
			
			delta_w=spectrum_match[:,1]-jacco_match[:,1]
			sigma_N=[]
			for w in range(len(delta_w)-4):
				w=w+2
				sigma_N.append(((delta_w[w+2]+delta_w[w+1]-delta_w[w-1]+delta_w[w-2])/2)**2)
			chi2_new=sum(sigma_N)
			if chi2_new<chi2_best:
				chi2_best=chi2_new
				best_parms = [FeH,input_ba,chi2_best]
	params_found.write(str(LEID)+','+str(temp_round)+','+str(grav_round)+','+str(best_parms[0])+','+str(best_parms[1])+','+str(best_parms[2])+'\n')
		
#pylab.plot(jacco_spec[:,0],jacco_spec[:,1],'k-')
#pylab.xlim([4505,4595])
#pylab.ylim([0.75,1.0])
#pylab.xlabel("Wavelength (angstroms)")
#pylab.ylabel("Normalized intensity")
#pylab.legend(('Observed spectrum','[Ba/Fe]=0.0','[Ba/Fe]=0.5','[Ba/Fe]=1.0'),'lower right')
#pylab.savefig("76041.pdf")
#pylab.show()
